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We present our first successful numerical results of 3D general relativistic simulations in which the 
Einstein equation as well as the hydrodynamic equations are fully solved. This paper is especially 
devoted to simulations of test problems such as spherical dust collapse, stability test of perturbed 
spherical stars, and preservation of (approximate) equilibrium states of rapidly rotating neutron 
star and/or corotating binary neutron stars. These test simulations confirm that simulations of 
coalescing binary neutron stars are feasible in a numerical relativity code. It is illustrated that using 
our numerical code, simulations of these problems, in particular those of corotating binary neutron 
stars, can be performed stably and fairly accurately for a couple of dynamical timescales. These 
numerical results indicate that our formulation for solving the Einstein field equation and hydro- 
dynamic equations are robust and make it possible to perform a realistic simulation of coalescing 
binary neutron stars for a long time from the innermost circular orbit up to formation of a black 
hole or neutron star. 
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I. INTRODUCTION 

The coalescence of binary neutron stars is one of the 
most promising sources for planned kilometer size laser 
interferometers such as LIGO |,|, VIRGO |], GEO |] 
and TAMA which will be in operation within the 
next five years. When a signal of gravitational waves 
is detected from binary neutron stars, it will provide not 
only the first chance to observe highly relativistic objects 
in a dynamical motion, but also a wide variety of physical 
information of binary neutron stars including their mass, 
spin, radius and innermost stable circular orbit (ISCO) 
. The signal of gravitational waves from such compact 
objects will be analyzed using matched filter techniques 
to extract the physical information. To apply this tech- 
nique, theoretical templates of gravitational waveforms 
are needed [Q. This fact has motivated an intense theo- 
retical effort for preparing such templates. 

The orbital evolution of binary neutron stars can be 
divided into three stages; inspiraling stage, intermediate 
stage, and dynamical stage. For the study of the inspi- 



raling stage in which their orbital radius is much larger 
than the stellar radius (R) and neutron stars are in quasi- 
periodic states gradually decreasing the orbital radius as 
a result of emission of gravitational waves, a post New- 
tonian approximation is powerful. Much effort has been 
paid for obtaining a template of high-order post Newto- 
nian corrections, providing many recent satisfactory re- 
sults §. 

When the orbital radius of the binary neutron stars de- 
creases to a few i? as a consequence of gravitational wave 
emission, the effect of the multipole moments of each neu- 
tron star induced by the tidal field from the companion 
star cannot be ignored. Even at this stage, the emission 
timescale of gravitational waves is still much longer than 
the orbital period. In this intermediate stage between 
the inspiraling and dynamical stages, the binary can be 
assumed to be in a quasi-hydrostatic equilibrium state. 
For the theoretical study, it is adequate to obtain the 
quasi-cquilibrium configuration taking into account the 
effects of the deformation of the neutron stars but as- 
suming that emission of gravitational waves is negligible. 
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In this research field, an effort has been paid recently, 
yielding gradually successful results 0-^. 

With the further emission of gravitational waves, the 
binary neutron stars approach the ISCO. Upon reaching 
the ISCO, they behave in a dynamical manner. For the 
theoretical study of such a dynamical stage, no approx- 
imation for general relativity is applicable because the 
system is dominated by highly general relativistic grav- 
ity and hydrodynamic effect. In this respect, numerical 
relativistic simulation is the only promising method for 
the theoretical study. 

Numerical relativity also plays an important role for 
a theoretical study on the origin of 7-ray bursts (GRBs) 
because the short rise times of the bursts imply that their 
central sources have to be relativistic objects @. Re- 
cently, at least some of GRBs have been turned out to 
be of cosmological origin ||ll|. In cosmological GRBs, the 
central sources must provide a large amount of the energy 
> 10^^ ergs in a very short timescalc of order msec-sec. 
It has been suggested that the merger of binary neutron 
stars could be a likely candidate for the powerful central 
source |0 . The typical scenario is based on the assump- 
tion that a system composed of a rotating black hole and 
a surrounding massive disk is formed after the merger. To 
clarify whether such scenario is correct, numerical simu- 
lations have been performed including effects of emission 
by gravitational radiation and/or neutrino and adopting 
a realistic equation of state for the neutron stars p2[ . 
So far, the results in the numerical simulations have not 
supported such scenario [Q; i.e., the evidence that the 
massive disk is formed around the black hole has not 
been found. However, all the simulations have been per- 
formed in the Newtonian or post Newtonian approxima- 
tions. Needless to say, general relativistic effects can play 
a very important role in the mergers between two neu- 
tron stars. To obtain the true answer, therefore, a fully 
general relativistic simulation is necessary. 

Much effort has been paid toward constructing a re- 
liable numerical relativity code which makes it possible 
to clarify the evolution of merging binary neutron stars 



and the gravitational waveform emitted by them. Sev- 
eral projects in the world such as those by Nakamura 
and Oohara and by Washington University group 

Ip^ are in progress, but no satisfactory results have been 
reported yet. 

To perform numerical simulation of coalescing binary 
neutron stars for a long time from the ISCO to forma- 
tion of a black hole or new neutron star, it is necessary to 
choose appropriate gauge conditions which make it pos- 
sible to perform the long-timescale simulation stably and 
to extract gravitational waves accurately. In a previous 
paper ]l6[ |, we performed fully general relativistic sim- 
ulations of coalescing binary clusters using coUisionless 
particles as the matter source of the Einstein equation. 
For the simulation, we used the approximate minimum 
distortion gauge and approximate maximal slicing condi- 
tions as the spatial and time coordinate conditions (see 
Sec II. B). We found that these gauge conditions are ro- 
bust enough to allow for stable and long-timescale sim- 
ulations of merging clusters as well as for the fairly ac- 
curate extraction of gravitational waves. In this paper, 
we perform simulations adopting the same gauge con- 
ditions and formulation for the Einstein equation and 
incorporating a solver of the relativistic hydrodynamic 
equations. We demonstrate the robustness of our for- 
mulation presenting successful numerical results of 3D 
hydrodynamic simulations. 

In particular, this paper is devoted to simulations of 
test problems. The purpose here is to ensure that numer- 
ical simulations of coalescing binary neutron stars for a 
long timescale from the ISCO to the formation of a black 
hole or neutron star are feasible in our formulation and 
numerical code. The test problems presented in this pa- 
per are a spherical dust collapse, stability test of spherical 
stars in equilibrium states, excitation of quadrupole oscil- 
lations of perturbed spherical stars, preservation of stable 
rapidly rotating stars in equilibrium states, and preserva- 
tion of corotating binary neutron stars in an approximate 
quasi-equilibrium orbit. We show that stable and fairly 
accurate simulations for these test problems are feasible 
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with our code, indicating the feasibihty of forthcoming 
reahstic simulations of coalescing binary neutron stars. 

The paper is organized as follows. In Sec. II, our 
formulation for solving the Einstein field equations as 
well as relativistic hydrodynamic equations are described. 
We also describe the gauge conditions and the numerical 
method employed in this paper briefly. In Sec. Ill, we de- 
scribe test problems which should be carried out to check 
the accuracy and performance of a numerical relativity 
code for solving the coupled equations composed of the 
Einstein and hydrodynamic equations. Then, we present 
the numerical results for the test simulations. In Sec. 
IV, we present numerical results of merger of corotating 
binary neutron stars as an example. Sec. V is devoted 
to a summary. Throughout this paper, we adopt the 
units G — c — Mq — 1 where G, c and Mq denote the 
gravitational constant, speed of light and the solar mass, 
respectively. Hence, the units of length, time, mass, and 
density are 1.477km, 4.927 x 10~*^sec, 1.989 x lO^^g, and 
6.173 X lO^^g/cm'^. Latin and Greek indices denote spa- 
tial components (1 — 3) and spacetime components (0—3), 
respectively. As spatial coordinates, we use the Carte- 



sian coordinates x'' ~ {x,y,z) with r = y^x^ + y'^ + . 
Sij{— 5'-^) denotes the Kronecker delta. 

II. FORMULATION 
A. Basic equations 

Our code solves the coupled equations of the Einstein 
equation and relativistic hydrodynamic equations. Our 
formulation for solving the Einstein equation has been 
described in detail in previous papers |17 lqjlq|. Since 



we adopt the same formulation only changing the matter 
source, we here only review the basic equations. 
We write the line element in the form 

ds^ = g^ydx^'dx'^ 

= {-a^ + PkP^)dt^ -t- 2(3,dx^dt + j.jdx'dx^, (2.1) 

where g^jy, a, /3* (Pi = ^ij(3^), and 7^ are the 4D met- 
ric, lapse function, shift vector, and 3D spatial metric. 



respectively. Following previous papers ||1^,|18|,|16[ , we de- 
fine the quantities to be solved in numerical computation 
as 



7 = det(7,,)^ei2*, 



-40, 



e '-7.y (i.e., det(7y) = 1), 



(2.2) 
(2.3) 
(2.4) 



where Kij is the extrinsic curvature, KjJ' its trace. We 
note that indices of Aij and/or are raised and lowered 
in terms of 7^ and 7*-' . In the numerical computation, 
we solve for 7^ , Aij , (f) and KjJ^ instead of 7^- and Kij . 
Hereafter, we use V^, Di and Di as the covariant deriva- 
tives with respect to g^^, jij and 7^, respectively. 

As the matter source of the Einstein equation, we 
adopt a perfect fluid. In this case, the energy momentum 
tensor is written as 



T^,i, ^ {p + pe + P)Uf,u^ + Pg^,y, 



(2.5) 



where p, e, P, and are the rest mass density, spe- 
cific internal energy density, pressure and four-velocity, 
respectively. Hereafter, we assume an equation of state 
in the form, P = {T — \)pe, where F is a constant. 

The hydrodynamic equations are composed of the con- 
tinuity, Euler and energy (or entropy) equations, which 
are derived from 



0, 



7/V. V = 0' 
u'^V.T/ = 0. 

We write their explicit forms as 

9tP* + 9.(p*w^) = 0, 
dt{p*Uk) + diip^iikv'') 



(2.6) 
(2.7) 
(2.8) 

(2.9) 



-ae 



.P 



whdkCt — UjdkP-' 



ae ^'^u^u^ _ - ■ 2ah{'w^ 
dkl ^ - - 



1), 



2wh 

(9te* -I- diie^.v'-) = 



w 



(2.10) 
(2.11) 



where = d/dx^, = pwe^'^, h = 1 + e + P/ p, w = 
au°, ilk = huk, — {peY^^we^'^ , and «*(= u^/vP) is 
written as 
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0:7*^ Uj 

whe^'f' 



(2.12) 
to 



In numerical simulation, we solve Eqs. (2.9)-(2.11 
evolve p*, Uk and e*. 

The volume integral of in the three hypersurface, 

EE J d^xp^, (2.13) 

denotes the total rest mass of the system, which should 
be conserved with time. Although the equation for has 
the same form as that for p*, the volume integral of is 
not conserved in the presence of shocks. This implies that 



using Eq. (2.11) we cannot obtain the correct solution 
when shock is formed. Thus, the artificial viscosity terms 



are added in Eqs. (2.10) and (2.11) for some problems in 
which shock is formed and plays an important role during 
the evolution of the system (see Appendix A). 

Once Ui is obtained, w{— au^) is determined from the 
normalization relation of the four- velocity, which can be 
written as 



= 1 



1 



rei 



^^-l 



(2.14) 



p*(we^ 

The Einstein equation is split into the constraint and 
evolution equations. The Hamiltonian and momentum 
constraint equations are written as 



D k 



where 



E = T^'^n^TLi, = p^hwe 



(2.15) 
(2.16) 



(2.17) 
(2.18) 



(— a,0), and Rij is the Ricci tensor with respect to 



Following our previous works, we write the evolution 
equations for the geometric variables in the form ]T8| , p^ 

= -2aA, + ^,kl3\, + ljkfi\, - (2.19) 



«( 



R 



+a(X/A,, - 2A,fcA/) 



(2.20) 
(2.21) 



(ai-/3'aO</'=^(-"^fe'+/5'fc), 
{dt-f3'di)K^'' =a[A,A' ^ 

-DkD^a + AT:a{E + S^^), (2.22) 

where Q i — diQ for an arbitrary variable Q, and 

5,, = T^'tm^.j = p.e-^'^iwhr^u.u, + e*%,P. (2.23) 



In calculating Rij and Rf!' in Eq. ( 2.20 ), we have terms 
of the type as S'^'^jikjj and 5^^^jk,ii- For evaluation of 
such terms, we define the auxiliary variable Fi = S^'-dijij 
Ip^ , ^ 16 1 and solve the evolution equation 



{dt - f3'di)F,= 2a{f^Ak,j + fijAk 



-2S^''aMA,,+S^'f3''jh,,,k 



(2.24) 



where hij = % - 6ij, and = 7*-' - 5'^. Then, 5''''%k,i] 



is evaluated as Fi 



We define the total angular momentum of the system 



as 



J = ^Ihn — <p [xAy, - yA.,,)e^^dS' 



xJy ~ yJx + ~ ^^x) 

87r " 



167r 
1 



Atjixdy -yda:)f- 



(2.25) 



where we use the Gauss's law and Eq. ( 2.16 ) to derive 
the final expression. We also use the following quantity 
to roughly estimate the angular momentum inside a co- 
ordinate radius r as 



J{r)= / d^xe''^{xJy - yJx). 

I <r 



(2.26) 



4 



B. Gauge conditions 

As in a previous paper |T^, we adopt an approximate 
maximal slice (AMS) condition and an approximate min- 
imum distortion (AMD) gauge condition as the time and 
spatial gauge conditions, respectively, for most of the 
simulations in this paper. In some test simulations car- 
ried out in Sec. Ill, we also use the zero shift vector 
gauge condition (3^ — for comparison between two re- 
sults obtained in different gauge conditions. 

To impose the AMS condition, we solve the following 
parabolic type equation for In a at each time step until 
an approximate convergence is achieved: 

dx Ina^ DkD'' Ina + {Dk hia){D'' Ina) - 47r(£; + S^^) 



(2.27) 



Here A denotes a control parameter and fa is a constant 
for which we assign a constant of 0(1). Assuming that 
the convergence is achieved and that the right-hand side 



of Eq. (2.27) becomes zero, the evolution equation for 
Kf}' can be written as 

[dt - = (2.28) 

Thus, if Kjj' is zero initially and the convergence is com- 
pletely achieved, the the maximal slice condition — 
is preserved. Even when the convergence is incomplete 
and KjJ' deviates from zero, the right-hand side of Eq. 



(2.28) enforces \K^\ to approach to zero in the local 

— 1/2 

dynamical timescale ~ . Hence, the condition 

Kj^ = is expected to be satisfied approximately. 

To impose the AMD gauge condition, we solve the fol- 
lowing simple elliptic type equations 



(2.29) 
(2.30) 



where A denotes the Laplacian in the flat 3D space, and 



= l&TTaJi + 2A,j{D^a ~ GaD^ct)) + -aDiK^^. (2.31) 

o 

From Pi and 7], we determine as 



P3 ^ ^ 



:P^ 



1 



(2.32) 



Namely, /3* satisfies an elliptic type equation in the form 
1 



(2.33) 

As we described in a previous paper [ [IGt , if an action 

/ = j d?x{da,j){dt%i)f''^^'. (2.34) 

is minimized with respect to we obtain the equation 
of a minimum distortion (MD) gauge condition for 

/?* as 

^jkD'b,l3^ + ^bjb,p' + Rjkf3'' = Sj, (2.35) 

where Rj^. is the Ricci tensor with respect to 7y . Thus, 
the equation for /3* in the AMD gauge condition is ob- 
tained by neglecting coupling terms between (3^ and hij 



in Eq. ( |2.35[ ). Since the neglected terms are expected 
to be small |jl^, we can expect that / is approximately 
minimized in the AMD gauge condition. 

The other benefit in the AMD gauge condition is that 
Fi is guaranteed to be small everywhere except in the 
strong field region just around a highly relativistic object 
JlGf . This implies that a transverse condition, S^^ di'^jk = 
0, approximately holds for 7.^ in the wave zone, helping 
the accurate extraction of gravitational waves near the 
outer boundaries of the computational domain. 

C. Initial value formalism 

Initial conditions are obtained by solving the con- 



straint equations ( 2.15 ) and ( 2.16 ). In this paper, we 
restrict our attention to initial conditions in which hij{= 
lij ^^ij) = ^^"i -^k' = 0. Then, the basic equations for 
obtaining the initial data are the same as those in [p8 16 
as described below. 

Using the conformal factor ip = e'^, Aij — tp^Aij and 
A*^ = ijj^A^^, the Hamiltonian and momentum constraint 
equations are rewritten in the form 

1 

8' 



(2.36) 



(2.37) 
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After we decompose Aij in the standard manner as 



= W,^j + Wj,, - -5,j5^'Wk,u (2.38) 

we set W, as 

W^z^^S, -i(x,. + Sfc,,x'=), (2.39) 
where x ^-nd Bi denote scalar and vector functions. Then, 



Eq. (2.37) can be decomposed into two simple elliptic 
type equations 



Ax = -SttJ.zV^- 



(2.40) 



Since Jiilj^{— p^iii) is non-zero only in the strong field 
region, the solution of the momentum constraint equation 
is accurately obtained. 

In addition to the constraint equations, we solve an 
elliptic type equation for a to impose the maximal slice 
condition initially. In the conformally flat 3D space, the 
equation is written in the form 

A(aV') = 2TTa^p^{E + 25/) + latfj-'^ A.jA'^ . (2.41) 

8 

From the initial condition, the total gravitational mass 
and angular momentum of the system at i = are cal- 
culated from 



1 



Jo = I d^x{xJy - yJx)i>^- 



An A 



(2.42) 
(2.43) 



Note that Eq. ( |2.25D reduces to Eq. ( ^.43D in the 3D 
space in which 7^ — Sij and KfJ^ = 0. 

D. Boundary conditions 

In this paper, we assume 7r-rotation symmetry around 
the z-axis as well as a plane symmetry with respect to 
the 2; = plane. Hence, we solve equations in a quadrant 
region L > x > —L and L > y, z > where L denotes 
the location of the outer boundaries. We impose the 
boundary condition in the y = plane such as 



Q{x, 0, z 
Q'^{x,0,z 
QaIx^O^z 
Q'ix,0,z 
Qz{x,0,z 
QABix,0,z 
Qaz{x,0,z 
Qzz{x,0,z 



= Q(-x,0,z), 
= -Q^(-a;,0,z), 
= ~Qa{-x,0,z), 
= Q^-x,0,z), 
= Qz{-x,0,z), 
= Qab{~x,0,z), 
= -QAzi-x,0, z), 

= Qzzi-X,0,z), 



(2.44) 

(2.45) 

(2.46) 
(2.47) 
(2.48) 
(2.49) 



where A and B — x or y, and Q, Q^{Qi), and Qij de- 
note arbitrary scalar, vector and tensor quantities, re- 
spectively. 

At the outer boundaries, we impose an approximate 
outgoing boundary condition for hij and Aij [|l8| as 



rQij {u) ~ const. 



(2.50) 



where we set u = at — e '*'r. (Even if we simply set 
u = t — r, the results do not change significantly.) More 



explicitly, Eq. ( 2.50| ) is rewritten in the form 
(5r ' 



1 



Qtj{t - 6t,r - 6r), 



(2.51) 



where St is a time step, and Sr — ae"'^'^St. Qij {t ~ St,r ~ 
Sr) is linearly interpolated from the nearby 8 grid points. 
We also note that the numerical results are not sensitive 
to the boundary condition of Aij and we have not found 
significant change in the numerical results even when we 
impose the boundary condition with the radial fall-off as 
0{r^^). A possible explanation for this result is that the 
spatial derivative of Aij , which appears in the evolution 
equations for Aij and Fi, does not play an important 
role for evolution of the system. On the other hand, 
numerical results and stability of the numerical system 
are significantly dependent of the boundary condition for 



jij. The condition defined by Eq. (2.51) is one of the best 
conditions among those we have tried so far. 

For (j) and KfJ^ , we impose the following boundary con- 
ditions at the outer boundaries 



Kk" - 0. 



(2.52) 
(2.53) 
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For Fi , we have tried a number of boundary conditions 
such as Fi = 0{r^'^), Fi = 0, and djFi =const at Xj = L, 
and have found that the results are weakly dependent 
of the boundary condition. We have found that the 
condition i^i = is preferable for a long-timescale nu- 
merical evolution, but for extracting gravitational waves 
near the outer boundaries, the condition Fi ~ 0{r^^) or 
djFi =const is preferable. In the spherical cases, we use 
the condition Fi = 0, but we change the condition case 
by case for other cases. 

It should be noted that all the outer boundary condi- 
tions described above are only approximate. This implies 
that numerical errors such as spurious back reflection and 
incoming of gravitational waves will be generated near 
the outer boundaries. For precise numerical simulation, 
apparently, we have to adopt more sophisticated outer 
boundary conditions as have been proposed and investi- 



gated by a couple of groups |21|. Imposing precise outer 
boundary conditions is one of important future issues. 



distribution changes in a dynamical timescale as in the 
case when the matter collapses to a singularity. Thus, we 
simply set 



St — min(C„ 



(2.56) 



32p, ' " ^' 

where C„i is a constant for which we choose 0.02 — 0.04. 
Note that the first term in the right-hand side of Eq. 
( ^.56 ) denotes the time for the collapse to a singularity 
of a spherical, homogeneous dust in the Newtonian limit. 
In the case when the density is so high that a black 



hole seems to be formed, the first term of Eq. ( 2.56 ) is 
smaller than the second term, but besides such a highly 
relativistic case, the second term is smaller and deter- 
mines the time step. We note that the hydrodynamic 
Courant condition is less severe than the geometric one 
p2|, so that we do not consider it. 



F. Brief summary of numerical methods and 
evolution scheme 



E. Grid and time step 

Throughout this paper, we use a uniform grid, i.e., 
5x ^ 5y = Sz ==const. We take (27V + 1, A^+l, iV+1) grid 
points in {x, y, z) direction, respectively (i.e., iV — L/ 5x). 

The time step 5t must satisfy the stability condition 
restricted by the Courant criterion for geometric vari- 
ables. If we neglect the other two directions, the geo- 
metric Courant condition in the direction is written 
as 



5t< [a7,^'/' + /3*]-ife. 



Since [a^ya ^^"^ + (3^] ^ is expected to be greater than unity 
for most cases, we simply set a geometric time step as 



5tg ~ CgSx, 



where Cg is a constant which we choose typically as 0.3. 
Also, St must be sufficiently small so that the matter 
distribution cannot change by a large fraction amount in 
one time step. The timescale is shortest when the matter 



For solving the equations for geometric variables as 
well as for determining the apparent horizon, the same 
methods as employed in previous papers [|l8|,^ are used 
here. The numerical method for solving the hydrody- 
namic equations is briefly discussed in Appendix A. 

The evolution scheme for the geometric and fluid vari- 
ables from a time slice at t to the next time slice at t + St 
is as follows (see the schematic figure 1): We put 7ij, Fi, 

(j), a, P'^, p*, iii, v^, and e* on t(°\ t^^\ ■ ■ ■ and Aij 
and K^'' on <(-i/2)^ t(i/2)^ ^(3/2)^. . i(«-i/2)^ ^^^^^ 

t^"^ denotes the coordinate time at the n-th time step, 
^(n+i/2) ^ (-^(n) ^ i("+i))/2, and n is a positive integer. 
Namely, we use the leapfrog method ||2^ for evolution of 
the geometric variables. 

For a given set of the geometric variables jij, Fi, 0, 
(2.55) a, (3^ and fluid variables p*, Ui, v\ and e* at t("\ and 
other geometric variables Aij and Kj!' at (the 
stage (1) of Fig. 1), first, Aij and K^}' are evolved from 



(2.54) 



^(n-i/2) i(n+i/2) ^gj^g ^j^g evolutiou equations ( |2.20| ) 
and ( p^ ) (the stage (2) of Fig. 1). Since the right- 
hand side of these evolution equations includes Aij and 
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Kf!^ which are defined only at i^" we use the hnear 
extrapolation method at each spatial grid point as 



I' 

3. r 



1 

1 , T 



to preserve the second order accuracy in time. We note 
that in adopting the linear extrapolation, we implicitly 
make use of the fact that the time steps at t^"^^) and 
t^"^ are approximately equal. 

Once we obtain (ii:^fe)("+i/2) and the ge- 

ometric variables 7ij, and Fi are evolved to (the 
stage (3) of Fig. 1). We also use the extrapolation 
method such as 

2 2 

because they are necessary to preserve the second order 
accuracy in time. 

Next, the hydrodynamic equations for p^, Uj, and 
are evolved to t*^"+^) . For solving the evolution equations, 
we use a Runge-Kutta method of second order accuracy 
psf ; i.e., the hydrodynamic equations are solved from 
to t("+i/2) in the first step (the stage (4) of Fig. 1), 
and using the fiuid variables defined both at t'"^ and at 
^(n+i/2) ^ those at are subsequently obtained in the 

second step (the stage (5) of Fig. 1). Since there appear 
7ij, (j), a, /3% w, h, and P in the right-hand side of the 
relativistic Euler equation ( 2.10D and we need those at 
^{n+i/2) jj-^ ^Yi^ second step for solving the equation, we 
use the extrapolation and interpolation, 

«("+l/2)^ 3^(n) _ 

2 2 



2 2 

P and h at are calculated from and d) 

at i("+i/2) obtained in the first step. (Note that 



and (Z?*)'-"''"^^ have not yet been obtained by this 
stage.) Once the fiuid variables are evolved to t("+^\ 



Eq. (2.14) is solved at each spatial grid point for ob- 
taining ■u;("+i) using the Newton-Raphson method [^ . 
Subsequently, we can obtain P and h at 

In the final step, we determine q;("+^) and (/3*)("+^) by 
imposing the gauge conditions (the stage (6) of Fig. 1). 
In solving their equations, we again use the extrapolation 
as 



(4)("+i) = 



3,7 

2 



(A,,)("+V2) 



1 . 7 

2 



(A,,)("-V2). 



Since all the quantities are evolved by this stage, we can 
derive at using Eq. ( ^.12| ) without any extrap- 

olation. 

III. TESTS AND RESULTS 

Our current priority in numerical relativity is to per- 
form simulations of the merger of binary neutron stars. 
Before carrying out such simulations successfully, it is 
necessary to confirm the accuracy and performance of 
our numerical code for many different problems. In par- 
ticular, the following issues have to be addressed: {i) the 
merger will take place for a couple of orbital periods from 
the time when the binary just enters inside the ISCO to 
formation of a black hole or neutron star. Can we carry 
out the simulation stably for a couple of orbital periods ? 
{ii) the final product of the merger will either be a black 
hole or a neutron star. If the merged object is unsta- 
ble against gravitational collapse, a black hole is formed. 
Can we judge the stability of the merged object against 
the gravitational collapse ? {in) the formation of a black 
hole will be signaled by the appearance of an apparent 
horizon. Can we determine the apparent horizon dur- 
ing the simulations ? (iv) can we extract waveforms of 
gravitational waves ? 

To answer these questions, we have performed simula- 
tions for a wide variety of test problems: 
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1. Spherical collapse of dust {P = 0) to a black hole: 
We compare the results with those obtained in a ID 
(spherical symmetric) simulation. We also check 
whether the apparent horizon can be found at a 
correct time and location. This test confirms that 
our code can simulate the formation of a black hole 
accurately. 

2. Stability of spherical stars: We check that the sta- 
bility of spherical stars can be judged in our code 
preparing both stable and unstable stars as initial 
conditions. We also check whether our code can 
provide a correct output on the fundamental radial 
oscillation frequency of perturbed spherical stars. 
This test is useful to confirm that we can determine 
the stability of the merged object against gravita- 
tional collapse. 

3. Quadrupolc oscillations of perturbed spherical stars 
and emission of gravitational waves: We give a 
quadrupole perturbation to a spherical stable star, 
and check whether we can obtain the frequency of 
the fundamental mode (/-mode) oscillation and ex- 
tract the waveform of gravitational waves near the 
outer boundaries. This test is useful to confirm 
that the extraction of gravitational waves near the 
outer boundaries is feasible. 

4. Preservation of rapidly rotating stars in equilibrium 
states: We check that (approximate) equilibrium 
states of rapidly rotating, stable stars can be pre- 
served for a couple of the rotation periods. The 
simulations are carried out choosing rigidly and 
rapidly rotating stars at mass-shedding limits. This 
test is useful to confirm that the coordinate twist- 
ing due to the rotation of the stars is sufficiently 
suppressed to a level adequate to carry out stable 
and accurate simulations in our AMD gauge condi- 
tion. 

5. Preservation of a corotating binary neutron star 
in a quasi-equilibrium state: We prepare a mildly 
relativistic corotating binary neutron star in an 



approximate quasi-equilibrium state obtained as- 
suming a conformally flat 3D geometry HJ]. Al- 
though we ignore hij which is necessary to ob- 
tain a true quasi-equilibrium configuration, it is at 
most a second post Newtonian quantity from the 
post Newtonian point of view [25| and for mildly 
relativistic binaries, the error is expected to be 
small. This test confirms that an approximate 
quasi-equilibrium state of a binary neutron star can 
be preserved for more than one orbital period. 

In the following subsections, we present the numerical 
results of the tests 1-5 separately. 

In the tests 2-5, we prepare the equilibrium and/or 
quasi-equilibrium states as the initial conditions adopting 
a polytropic equation of state P = Kp^ with F = 5/3 
and/or 2. For F = 5/3 and 2, we fix K as 10 and 200/7r, 
respectively, to mimic a relation between the rest mass 
Af* and the central density pc for neutron stars. 

In Fig. 2, we show and the circumferential ra- 
dius R of the spherical equilibrium stars as a function 
of Pc for F = 5/3 and 2, respectively. For F — 5/3, 
Af* {Mg) reaches a maximum value ~ 1.552 (1.487) at 
Pc =i 1.85x10^3^ and for F = 2, Af* {Mg) ~ 1.435 (1.306) 
at Pc — 5.0 X 10^"'. Beyond the critical densities, the star 
is unstable. 

Although we adopt particular units fixing K , the mass, 
length, and density may be rescaled using the following 
rule 

M^{Mg) M^C'/^iMgC"'^), R -> i?C"/^ 

Pc ^ PcC"" and J ^ JC" for if ^ if C, (3.1) 

where n = 1/(F — 1) and C is an arbitrary constant. 
Namely, the invariant quantities are only the dimension- 
less quantities such as M^R-^^^ {MgK^'^/^), RR-'^I'^, 
PcK"", M,/R{Mg/R), and J/Ai^ (J/M^). 

A. Test 1: Spherical collapse of dust 

We consider a time symmetric, conformally flat initial 
condition for the dust sphere, and adopt the following 
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density profile 

p* ^ A 



1 + cxp 



,2 s 1-1 



(3.2) 



where we choose tq — 4Mg and 6r'^ = O.lSMg. A is 
adjusted so that the gravitational mass of the system is 1 
{A ~ 4.287 X 10^'^). In this test, we assign the negligible 
specific internal energy and pressure, i.e., e <^ 1 and 
P p. Throughout this subsection, every quantity is 
shown in the unit Mg — 1 (and G = \ = c). 

We perform the 3D numerical simulation as well as 
ID (spherical symmetric) simulation and compare two 
results. The simulations are carried out in two different 
spatial gauge conditions, the MD and /?' = gauge con- 
ditions. We note that in a spherical symmetric case with 
a conformally flat initial condition, the AMD gauge con- 
dition is identical with the MD gauge condition because 
the condition dt^fij = holds in both cases throughout 
the whole evolution. Hence, the comparison can be done 
without any coordinate transformation. We use the max- 
imal slice condition in the ID case and AMS condition in 
the 3D case. Since Kj!' can be kept nearly equal to zero 
in the AMS condition, we can consider them as the same 
slice conditions. We have used a grid with iV = 60 and 
5x = 0.15 in the 3D simulation. 

First, we show numerical results for the (3^ = case. In 
Fig. 3, we show a and at r = as a function of time. 
The dotted and solid lines are the ID and 3D results, 
respectively. In Fig. 4, we also plot a along the x-axis at 
selected times {t = 12.0, 17.2 and 20.4) for the ID (the 
solid lines) and 3D (the filled circles) results, respectively. 
We note that in the ID simulation, the apparent horizon 
is found a,t t ^ 18.6. In the 3D simulation, it is also 
found nearly at the equal time and the same location. 
These results show that the ID and 3D results agree well 
at least up to the formation of a black hole and confirm 
that it is possible to investigate a black hole formation 
in a spherical symmetric spacetime accurately with our 
formulation and numerical code. 

In the late phase {t > 23), the accuracy of the 3D 
results deteriorates. This is because the error increases 



in 7y rapidly as a consequence of the so-called horizon 
stretching around the location of the apparent horizon. 

Next, we present numerical results in the MD (and 
AMD) gauge conditions. In Fig. 5, we show a and p, 
at r = as a function of time. In this case, the ID (the 
dotted lines) and 3D (the solid lines) results agree well 
by t ~ 20. Fortunately, we could determine the apparent 
horizon because it is formed t ^ 18.6, but the accuracy 
deteriorates soon after the formation. This illustrates 
that the MD gauge condition is less appropriate than 
P'' ~ gauge condition for evolution of the late phase of 
the gravitational collapse in the simulations performed 
under identical grid number and spacing. The reason 
is apparently related to the drawback of the MD gauge 
conditions pointed out in : In the MD or AMD gauge 
conditions with K)} ~ 0, the coordinates spread outward 
and the proper distance between two neighboring grids 
increases, i.e., > (/?'" > 0) and dtcj) > (cf. Eq. 
( ^.2l| )), around r = during gravitational collapses. As 
a result, the black hole forming region cannot be well 
resolved. In Fig. 6, we show the time evolution of at 
r = 0. It changes from 0.17 to 0.9 by t = 20. Thus, the 
proper distance between two neighboring grids increases 
by a factor e2("#'(*=2O)-0(t=o)) _ 4_ ^^le f3' = gauge 
condition with KfJ' = 0, on the other hand, (/'(a;*) is 
constant in time and hij for r is nearly equal to 0, so 
that the physical grid spacing in the MD gauge conditions 
is 4 times as large as that in the = gauge condition 
around r = a.t t ^ 20. 

To overcome the coordinate spreading effect without 
changing the gauge condition, we have to take a large 
number of grid points around the black hole forming re- 
gion. The solution to this problem may be to adopt an 
adaptive mesh refinement technique p6[ | in which we can 
improve the resolution around the black hole forming re- 
gion effectively. However, since the spreading factor in- 
creases rapidly as shown in Fig. 6, we have to improve 
the resolution also quickly when using such a technique. 
We do not think it a good idea to simply rely on such a 
technique. We consider it necessary to modify the gauge 
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condition. A strategy of the modification and the numer- 
ical experiment will be presented in one of forthcoming 
papers [ p7| . 

B. Test 2: Stability of spherical stars 

In this second test, we prepare spherical equilibrium 
stars of a polytropic equation of state of F = 5/3. We 
use two stars for the test. One is a stable star of pc = 
10-3, Af, = 1.499, and Mg = 1.440, and the other is 
an unstable star of pc = 2.4 x 10^"^, = 1.542, and 
Mg — 1.478. In both cases, the total rest mass is slightly 
less than the maximum value. In Fig. 2, we show with 
the filled circles and the open circle the locations of the 
two stars and the star at the critical density of stability, 
respectively. 

Here, we again use two spatial gauge conditions: AMD 
gauge and = gauge conditions. In Fig. 7, we show 
the time evolution of the density p(= p^.e~^'^/w) and a 
at r = as a function of time for a stable configuration 
of Pc = 10^3. These simulations were performed with 
= 50 and Sx = 0.4. The solid and dotted lines are 
the results in the /?* = gauge condition, and the dashed 
line is in the AMD gauge condition. The dotted line also 
denotes the result for the case in which we give a per- 
turbation to the equilibrium configuration by reducing 
K of 0.5% initially. In other cases, we give the equilib- 
rium state without any changes. These figures clearly 
illustrate the feasibility of our code to judge the stability 
of the spherical stable star and to preserve it stably at 
least in a few oscillation periods irrespective of the gauge 
conditions. 

As shown with the approximate perturbation analysis 
in appendix B [p8| , the period of the hmdamental radial 

— 1/2 

oscillation is ^ 10.5pc . For the pressure depleted case 
(the dotted line), the oscillation period is clearly recog- 
nized in Fig. 7. This shows that the frequency of the fun- 
damental radial oscillation can be computed accurately 
in our code. 

In Fig. 8, we show the time evolution of p and a at r = 



for an unstable configuration of pc = 2.4 x 10"^. The 
simulations are performed with A^ = 50 and 6x = 0.3. 
We prepare pressure depleted initial conditions in which 
we reduce K by 0.5% and 0.1%. As in the stable case, 
we use both the /3' = and AMD gauge conditions, and 
the numerical results are denoted by the solid and dashed 
lines, respectively. As shown in Fig. 8, the unstable star 
collapses into a black hole. As one would expect, the star 
collapses more quickly in the case when the depletion 
factor of the pressure is larger. These results confirm 
that our code provides correct results. 

The time evolution of p and a at r = in the two 
spatial gauge conditions approximately agrees with each 
other except for the late phase of the gravitational col- 
lapse during which the coordinate spreading effect is se- 
vere in the AMD (and MD) gauge conditions The 
overall agreement (except for the late phase) is an im- 
portant test of consistency because p and a at r = are 
both gauge independent quantities. 



C. Test 3: Quadrupole oscillations of a perturbed 
spherical star 



For the third test, we prepare stable spherical stars in 
equilibrium states as in test 2. We adopt F = 5/3 and 2 
in this test, and choose the stars in which pc — 10^^ for 
F = 5/3 and pc = 3 x 10^^ for F = 2, respectively. For 
F = 2 and Pc = 3 X 10"^, the mass and circumferential 
radius are 1.25 and 6.99, respectively (see filled circles in 
Fig. 2). We perform the simulations with A^ = 50 and 
Sx = 0.5 for F = 5/3 and with A^ = 50 and 5x = 0.25 
for F = 2, respectively. The test is performed in the 
AMD gauge condition and at the outer boundaries, we 
set F, = 0(r-3). 

As the source of the initial quadrupole oscillation, we 
give a velocity perturbation of type 



u^{t = 0) 
or of "x" type 



A 




u.^{t = 0) = A\lj^{~y,~x,0), 



(3.3) 



(3.4) 
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where A is a constant which is set to be 0.082. 

As a result of a time- varying mass quadrupole moment 
of the star, gravitational waves are emitted. This allows 
to check the feasibility of the gravitational wave extrac- 
tion near the outer boundaries. To extract gravitational 
waveforms, we define |16| 



For r = 5/3 and pc = lO'^ (i.e., Mg ~ 1.44 and 



h+ = r{'y^^ ~ %a)/2. 



(3.5) 
(3.6) 



along the z-axis. Since we adopt the AMD gauge condi- 
tion and prepare initial conditions in which Fi — 0, hij is 
approximately transverse and traceless in the wave zone 
p6[ . As a result, and are expected to be appropri- 
ate measures of gravitational waves. They are also useful 
to find the maximum amplitude of gravitational waves 
because we treat the problems in which the amplitude is 
maximum along the z-axis. 

In Figs. 9-11, we show a root mean square radius 

1/2 



(3.7) 



as a function of time for the + and x mode perturbations 
for F — 5/3 (Figs. 9 and 10), and for + mode pertur- 
bation for F = 2 (Fig. 11), respectively. In the case 
of the + mode oscillation, Xims (the solid line) and t/r,ns 
(the dotted line) oscillate with a characteristic period. 
For F = 5/3 and 2, the oscillation period is ^ 5.5 pc 

— 1/2 

and ^ 4.7pc , respectively. As shown by perturbative 
studies on stellar pulsations p0| , the angular frequen- 
cies of the /-mode oscillation for F = 5/3 and 2, and 
Mg/R ~ 0.1 are approximately written as 



1.44a 



1.22 




for r = 5/3. 



for F 



(3.8) 



Or, according to an empirical formula pi| , the angular 
frequency for the stars of Mg/R ^ 0.1 is also written, 
irrespective of F, as 



u) ~ 0.012 + 0.93 




(3.9) 



R ~ 13.1), the oscillation period is ~ 5.5p, 



-1/2 



both 



formulas. For F = 2 and pc = 3 x 10"^ (i.e., Mg ~ 1.25 
and R = 6.99), the period is ~ 4.7 pc^^'^ from Eq. (jsj) 



— 1/2 A ^ 

and ~ 5.0pc from Eq. (3.9). Thus, our numerical 



results agree with these perturbative results fairly accu- 
rately. 

On the other hand, the oscillation periods of Zj-ms (the 
dashed line) for the + mode perturbation and of all the 
components of xl^-^^ for the x mode perturbation are 



roughly - lOpc for F = 5/3 and - 7pc '''' for F = 2. 
These are in agreement with the periods of the funda- 
mental radial oscillation and different from those of the 
/-mode oscillation. This reflects the fact that they are 
not relevant for the /-mode oscillation at linear order. 

In Figs. 12 and 13, we show and /ix at Zobs — 24.5 
(the solid lines) and 19.5 (the dashed lines) for F = 5/3 
as a function of the retarded time for the + and x mode 
perturbations. In Fig. 14, we also show and ft-x at 
^obs — 12.25 (the solid lines) and 9.75 (the dashed lines) 
for F = 2 for the + mode perturbation. As the oscilla- 
tion frequency of the /-mode shows, the wavelength of 
gravitational waves is several times larger than Zobs, so 
that we expect the extracted waveforms to be different 
from the asymptotic waveforms. However, they appear to 
constitute a fair description of the asymptotic waveforms, 
at worst qualitatively, because (a) the solid and dashed 
lines approximately agree, i.e., '^^x — lyy and 'yxy behave 
approximately as /(t — z)/ z along the z-axis where / de- 
notes a generic function, ( h) the oscillation period agrees 
well with that of the /-mode, ( c) the order of magnitude 
of and h-^ roughly agrees with that derived by the 
quadrupole formula, i.e., Mgv"^ ^ MgA^/R ^ lO^'', and 



-1/2 



(d) for the perturbation given by Eq. (3.3), h X remains 



nearly zero, and for Eq. (3.4), remains nearly zero, 
as expected by the quadrupole formula. These facts sug- 
gest that and hx represent (at least approximately) 
gravitational waves emitted by the stellar oscillation and 
indicate the ability of our code to extract gravitational 
waves near the outer boundary. We expect that a simula- 
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tion of higher resolution would produce even more precise 
gravitational waveforms. 

D. Test 4: Stability of rapidly rotating stars 

To carry out this test, we prepare rapidly and rigidly 
rotating neutron stars in (approximate) equilibrium 
states. We adopt F — 5/3 and 2, and select stable com- 
pact neutron stars. 

To prepare the rotating stars in approximate equilib- 
rium states, we use a conformal flatness approximation. 
Then, the geometric and hydrostatic equations for solu- 
tions of equilibrium states are described as 

o 

A{ai;) = 2TTatP^[ph{3w'^ - 2) + 5P] 

+ ^P''P'L,,Lki, (3.11) 



r = 2. As the density of the rotating stars increases (i.e., 
Pt), the coincidence becomes worse gradually, but 



Pc 



> 



~ IGnaphwUj, 



ah 

— = const. 

w 



where 



h = l + KTp^-^/{T-l), 

ia\ 6 J 



(3.12) 
(3.13) 

(3.14) 
(3.15) 
(3.16) 
(3.17) 



and r2 denotes the angular velocity of the rotation. 

Although the solutions obtained from Eqs. ( |3.1C| )- 
( |3.13 ) are not exact, we can still expect that they are 
excellent approximate solutions as illustrated in [p^ . As 
shown in Fig. 15, this is the case: We show the grav- 
itational mass Mg as a function of the central density 
Pc for rotating stars at mass-shedding limits obtained 
from exact equations (the solid lines) |Q and obtained 
by the conformal flatness approximation (the circles). We 
have found that the sequences of the circles almost co- 
incide with the solid lines for mildly relativistic stars of 
Pc < Pt where pt ^ 0.0015 for V = 5/3 and pt - 0.003 for 



even for pc ^ 2pt at which the star seems to be unstable, 
the difference between the sequences of circles and solid 
lines is small. 

We perform the simulations for the rotating stars 
marked with the filled circles in Fig. 15. The relevant 
quantities for the rotating stars are shown in Table I. We 
adopt Sx = 0.434 for F = 5/3 and 5x = 0.202 for F = 2 
as the grid spacing. With these grid spacings, the major 
and minor axes of the stars are initially covered by 40 
and 23 — 24 grid points, respectively. 

To suppress the coordinate twisting, the AMD gauge 
condition is adopted. The shift vector (3'^ determined 
in this gauge condition for the conformal 3D metric 
{jij = Sij) agrees with that obtained from Eq. (3.12). 



This implies that when a simulation is started, the gauge 
condition is identical with that used for obtaining the ap- 
proximate equilibrium states. Therefore, we can use the 
approximate equilibrium states as the initial conditions 
without any coordinate transformation. 

In this test, Fi is set to be zero at the outer boundaries. 
The grid number N is set to be 54 typically. We also 
adopted N — 76 without changing the grid spacing in 
some of the following simulations, but we did not find 
significant difference in the results. This indicates that 
the outer boundary condition adopted here is adequate. 

We have considered two kinds of initial conditions: In 
one case, we use the approximate equilibrium configu- 
rations without any change. In the other case, we ini- 
tially decrease the pressure by reducing K of 1% (i.e., 
AK/K = 1%). Even in the pressure depleted case, the 
stars should be stable because the gravitational mass is 

3 — 5% smaller than the maximum mass of the rotat- 
ing stars (see Fig. 15). In Fig. 16, we show snapshots 
of the density contours lines for p^ and the velocity field 
for {v^, v^) in the equatorial plane (left) and in the y = 
plane (right) at selected times for AK = and for F = 2 
as an example (for F — 5/3, we have found that similar 
figures can be drawn). We also show a and p at r = as 
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a function of t/P in Fig. 17, Xrms and Zims as a function 
of t/P in Fig. 18, and J/ Jo as a function of i/P in Fig. 
19 for AK = (the solid lines) and AK/K = 1% (the 
dotted lines) and for F = 5/3 and 2. Because of numer- 
ical dissipation at the stellar surface, the total angular 



momentum of the system (J in Eq. ( 2.25 )) decreases by 
^ 2% by i ^ 2P in all the simulations (see Fig. 19). As 
a result, a at r = and xj^g decrease and p at r = 
increases with the time evolution. Also, the stars suffer 
slight non-axisymmetric (quadrangular shape) deforma- 
tion near the surface (see Fig. 16) because we use the 
Cartesian coordinates for rotating stars of a spheroidal 
shape. However, besides these slight changes, the stars 
remain almost in the stationary states for more than two 
rotational periods (i.e., by the time when we terminated 
these simulations). 

Since the initial conditions adopted are only approx- 
imate equilibrium states, some oscillations are induced 
with time evolution. However, the amplitude is very 
small and cannot be distinguished from the numerical 
errors. Thus, the states of the rotating stars are close to 
the true equilibrium ones. Actually, the absolute value of 
each component of hij remains small and of order ~ 0.05. 
These results re-confirm that the conformal flatness ap- 
proximation is really a good approximation for obtaining 
axisymmetric rotating stars in equilibrium states. 

For AK/K ~ 1%, a and p at r = and xl^^ oscillate 
with the time evolution. The oscillation period is roughly 
0.9P for both F = 5/3 and 2. We deduce that the period 
denotes a fundamental quasi-radial oscillation period of 
the rotating stars. 

We emphasize that the simulations can be stably car- 
ried out for more than two orbital periods without any 
instabilities and with hij being kept small. These re- 
sults clearly demonstrate that our formulation for solving 
the Einstein equation is robust even for systems of non- 
zero angular momentum and that the coordinate twist- 
ing is sufficiently suppressed to a level adequate for long- 
timescale simulations in our AMD gauge condition. 

In this paper, we only perform test simulations us- 



ing stable rotating stars. It is very interesting to per- 
form simulations adopting rapidly rotating supramassive 
neutron stars to investigate the stability and the fate of 
unstable stars |^8| . Such rapidly rotating supramassive 
neutron stars may be frequently formed as a result of 
accretion onto a neutron star from a companion star in 
normal binary systems p9| . However, such simulations 
are beyond the scope of this paper. More detailed anal- 
ysis of the stability of rapidly rotating stars and of the 
final fate of unstable stars will be presented in p^]. 



E. Test 5: Corotating binary in an approximate 
quasi-equilibrium state 



In the final test, we adopt a mildly relativistic coro- 
tating binary neutron star in an approximate quasi- 
equilibrium state as initial condition. As mentioned in 
Sec. I, our purpose in future is to carry out simula- 
tions of coalescing binary neutron stars from the ISCO 
to formation of a black hole or new neutron star. In the 
early stage in which the binary is near the ISCO, the 
radial velocity of each star is expected to be small, and 
we can consider the binaries as on approximate quasi- 
equilibrium orbits rather than on plunging orbits. As 
the radial velocity gradually increases, the orbit changes 
from the near inspiraling one to the plunging one. For a 
successful simulation of a coalescing binary neutron star, 
therefore, the nearly quasi-equilibrium orbit has to be 
maintained at least for 1 orbital period stably. In this 
test, we show this feasible with our code. 

The configuration of a binary neutron star is again 
obtained in the assumption of the conformally flat 3D 
metric and the maximal slice condition Kj}' = [p4[ . 
We prepare a binary in which the surfaces of two stars 
come into contact. We adopt F = 5/3 in this test. The 
equation of state for such small F < 2 is not so stiff that 
binary neutron stars in contact orbits are stable against a 
hydrodynamic instability [p4[ . Since the binary neutron 
star is not very compact (see Table II) and is also far from 
the ISCO, it can remain on the stable orbit for a time 
comparable with the emission timescale of gravitational 
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waves. 

Approximate quasi-equilibrium states of binary neu- 



tron stars are obtained solving Eqs. ( |3.10| - |3.13D . The 
equations are solved using a numerical method similar to 
that employed in |Q. 

In Fig. 2, we denote, by the cross (of higher density, 
Pmax — 10^'^), the relation between the maximum den- 
sity Pmax and half of the total rest mass M*/2 for the 
corotating binary neutron star in an approximate quasi- 
equilibrium state which is used as an initial condition 
of the simulation. We compared the numerical results 
on the relation between the maximum density and rest 
mass with those by Baumgarte et al. |2^, and found that 
they agree within 2% error. In Table II, we also list the 
relevant quantities of the binary neutron star. We define 
an orbital radius a = Al^^fl''^^^ , and we define an ap- 
proximate ratio of the emission timescale of gravitational 
waves to the orbital period as 



Rr 



1287r V M, 



5/2 



1287r 



(Mgn) 



-5/3 



(3.18) 



where we have used the Newtonian expression of the en- 
ergy and orbital period, and the quadrupole formula for 
the energy luminosity of gravitational waves Q . Since 
our purpose is to check the feasibility of our code to pre- 
serve the approximate quasi-equilibrium state stably, this 
test should be performed for a binary in which Rr > 1. 

In the numerical simulation, we adopt the AMD gauge 
condition to sufficiently suppress the coordinate twisting. 
As we discussed in Sec. III. D, the gauge condition at < = 
is identical with that used for obtaining the approximate 
quasi-equilibrium state, so that we do not have to carry 
out coordinate transformation at t = 0. 

In this subsection, we adopt the grid spacing as Sx = 
0.927, varying the grid number N as 76 and 116. With 
this setup, the major diameter of each star is covered by 
30 grid points initially. For = 76 and 116, i(= N6x) 
is only ^ 0.2 and 0.3 times as large as the wavelength of 
gravitational waves Agw = vr/f2. For a precise simulation, 
the outer boundary should be located in the wave zone 
so that L ^ Agw. This is because we should impose an 



outgoing boundary condition for 7^ and Aij , and in the 
system of binary neutron stars, the existence of gravi- 
tational waves plays an important role for evolution of 
the system. We have found that the present incomplete 
treatment for the outgoing boundary conditions seems 
to produce numerical errors (see below further discus- 
sion). However, for imposing the boundary condition in 
the wave zone with the uniform grid, it is necessary to 
take a very large grid number A^ > 400. That is not feasi- 
ble in the present computational resources on supercom- 
puters. In this paper, we perform simulations overlooking 
the deficits. Thus, the radiation reaction of gravitational 
waves is not precisely taken into account in the follow- 
ing simulation. A large simulation of L > Agw is one of 
future issues. 

In Fig. 20, we show snapshots of the density contour 
lines for and velocity field for (u^,w^) in the equato- 
rial plane at selected times for A^ = 116 while in Fig. 21, 
we show a:;*jjjg as a function of t/P. We note that if the 
binary remains in a quasi-equilibrium state, the orbital 
period should be kept to ~P, the curves for x^ns and 
Urms should be complete sine curves, and Zj-ms should be 
a constant. From Figs. 20 and 21, it is evident that the 
state of binary neutron star fluctuates from the initial 
state with time. The fluctuations seem to be mainly due 
to the numerical error discussed below. Besides the spu- 
rious numerical effect, the binary neutron star is kept in 
an approximate quasi-equilibrium circular orbit for ^ 1 
orbital periods stably. 

There appear to be two candidates for the numerical 
error. One is the numerical dissipation of the angular 
momentum at the stellar surface p5[ | which was already 
pointed out in test 4. This is a simple consequence of 
insufficient resolution. As a result, the orbital radius de- 
creases. The other candidate is the incomplete treatment 
of the outgoing boundary conditions. As mentioned be- 
fore, we impose an approximate outgoing boundary con- 
dition for and Aij not in the wave zone. From this 
incomplete treatment, the angular momentum seems to 
go out and come in inaccurately from the outer bound- 
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aries. As a result, the orbital radius increases and de- 
creases. To illustrate the fact, J/ Jo as a function of t/P 
is depicted in Fig. 22. From this figure, we can recog- 
nize that the angular momentum increases and decreases 
with time by ~ ±3% of the total angular momentum. 
We note that the angular momentum should be mono- 
tonically dissipated by gravitational radiation by about 
2% of the total angular momentum in one orbital period 
according to the estimation by the quadrupole formula. 
However, the numerical results do not reflect this effect 
accurately. 

Despite these errors, however, the binary neutron star 
is preserved in the approximate quasi-equilibrium orbit 
for more than one orbital periods stably. This seems to 
imply that we are choosing an adequate gauge condition 
and formulation of the Einstein equation. We therefore 
expect that if we could perform a simulation taking a suf- 
ficient number of grid points to impose the outer bound- 
ary conditions in the wave zone and to improve the res- 
olution, it would be possible to perform the simulation 
not only stably, but also accurately. 

During the evolution, hij deviates from zero and grad- 
ually reaches a finite amplitude. The maximum absolute 
value of each component is of order ^ 0.05 by ^P. (This 
is roughly consistent with what is inferred from a post 
Newtonian study |p5| in which hij has magnitude of or- 
der ~ [Mg/aY-) Since hij becomes non-zero, it could 
slightly affect the quasi-equilibrium configuration of bi- 
nary neutron stars. However, this does not seem to be 
a serious perturbation for the present mildly relativistic 
case. 

In Fig. 23, we show /i+ and /ix as a function of t/V . 
The solid and dotted lines are extracted at Zobs — 106.6 
and 85.3 for N = 116, and the dashed line at Zobs — 69.5 
for N = 76. As mentioned before, we extract them at 

0.2 — 0.3Agw, so that they only indicate approximate 
asymptotic waveforms. Nevertheless, we find that they 
behave in the periodic manner and that the period ap- 
proximately agrees with the orbital period. Furthermore, 
the solid and dotted lines agree very well, implying that 



the waves propagate at the speed of light. These facts in- 
dicate that the approximate waveforms constitute a fair 
description of the asymptotic ones. 

There are, however, differences from what is expected 
in the asymptotic waveforms. The first one is found in the 
long wavelength modulation (apparent especially in 
which should not appear in the correct waveforms. The 
modulation is larger for the simulation with smaller iV, 
indicating that it is caused in the outer boundaries. The 
second difference is found in the irregular, non-periodic 
waveforms in the early stage for t < Zobs • These waves are 
emitted because we set an approximate quasi-equilibrium 
state as the initial condition neglecting hij . The first ra- 
diation will be the relaxation of the system to the true 
quasi-equilibrium state. To avoid this shortcoming, we 
have to prepare more realistic quasi-equilibrium states 
adopting a formalism in which hij is appropriately taken 
into account. The third difference is found in the am- 
plitude: For the case where two point particles of equal 
mass are in a circular orbit, the maximum amplitudes of 
/i+ and ft-x hi the post Newtonian approximation can be 
written as IM 
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(3.19) 



where x = (MgD,)"^^^ . Using this formula, the maximum 
value of h+ and hx in the binary of Mg — 2.92 and 
X = 0.0887 is 0.24. (Note that in the quadrupole formula, 
it is 0.26.) Since the convergence of the post Newtonian 
formula for x 0.1 is not very good, we should take 
into account an error of ~ 10%. Even if wc consider such 
error, the amplitude obtained in the numerical simulation 
is found to be somewhat smaller than that derived from 



Eq. ( p.l9D . The reason seems again due to the incomplete 
treatment of the outgoing boundary condition at L ~ 
0.2-0.3Agw. 

In summary, we can perform simulations of corotating 
binary neutron stars in an approximate quasi-equilibrium 
state stably and fairly accurately, and extract gravita- 
tional waves with ^ 10% error even in the present re- 
stricted grid numbers. However, to improve the accuracy 
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for gravitational waveforms, it is necessary to adopt more 
sophisticated bomidary conditions to take a larger 
number of grid points and to prepare a more realistic 
quasi-equilibrium state as the initial condition. These 
issues will be addressed in future simulations. 

IV. MERGER OF COROTATING BINARY 
NEUTRON STARS: A DEMONSTRATION 

In this section, we show numerical results of merger be- 
tween two neutron stars as demonstration that such sim- 
ulations are feasible. We again adopt corotating binary 
neutron stars of F = 5/3 in contact and in approximate 
quasi-equilibrium orbits. Two binaries shown in Table II 
(or in Fig. 2 with crosses) are prepared as the initial con- 
ditions: One is the same as that adopted in Sec. III.E, 
and the other is less relativistic one in which the initial 
maximum density of each star is 6 x lO^**. To accelerate 
the merger, we artificially reduce the angular momentum 
by 5% in the initial stage; i.e., we reduce the initial values 
of Ux and Uy by 5% uniformly. In this section, we adopt 
N = 116. 

In Figs. 24 and 25, we show snapshots of the density 
contours lines for and the velocity field for (w^,?;^) 
in the equatorial plane for Pmax(i = 0) = 10^"^ and 
6 X 10"'', respectively. Since we artificially reduce the 
angular momentum of the system at t = 0, the neutron 
stars approach each other to merge soon after the simu- 
lations start. For both types of initial data at t > 0.5P 
(P denotes the orbital period of the quasi-equilibrium 
states without the angular momentum depletion), the 
neutron stars begin to merge forming spiral arms, and 
by t ~ 1.5P, new neutron stars which are nearly axisym- 
metric are formed. In Figs. 26 and 29, we show snapshots 
of the density contour lines for in the j/ = plane at 
t = 1.62P for p,nax(t = 0) = 10-^, and at t = 1.59P 
for Pmax(^ = 0) = 6 X 10"^. In Fig. 28, we also show 
the angular velocity il = {xv^ — yu^)/(x^ -I- y'^) along 
the X and y-axes in the equatorial plane at t = 1.62P for 
Pmax(t = 0) = 10-3 and at t = 1.59P for p^^^{t - 0) = 
6 X lO-**, respectively. These results show that the final 



products are rapidly and differentially rotating, highly 
flattened neutron stars. For the case Pmax{t = 0) = 10^^, 
the central density of the merged object for t ~ 1.6P is 
^ 1.4x 10~3, which is nearly the maximum allowed den- 
sity along the sequence of stable neutron stars of K — 10 
and r — 5/3 (cf. Figs. 2 and 15). This indicates that the 
new neutron star is expected to be located near the crit- 
ical point of stability against gravitational collapse and 
that a black hole could be formed in the merger of more 
massive neutron stars. 

In Fig. 29, we show fraction of the rest mass inside a 
coordinate radius defined as 

as a function of time for Pmax(^ = 0) = 10"^ (the solid 
lines) and pmax(i = 0) = 6 x 10"'' (the dashed hues). 
The figures show that ~ 96% and ~ 97% of the matter 
are inside r = 18(~ 6Afg) for pmax(i = 0) 10"'^ and 
r = 20.4(~ 7.5Afg) for pmax(i 0) = 6 x 10"* forming 
the new neutron stars. Thus, only a small fraction of the 
matter can spread outward to form a halo around the 
central objects. We also find that J(r = 18) ~ 0.85Jo 
for pmax(t = 0) = lO"'"^ and J(r = 24.4 9Mg) - 0.85Jo 
for Pmax(^ = 0) = 6 X 10"''. In both cases, J/M'^ of the 
newly formed neutron stars seems roughly ~ 1 , and that 
appears to be one of the reasons for which the neutron 
stars do not collapse to be a black hole. 

In Fig. 30, we show /i+ and as a function of 
t/P for the two cases. For Pmax(^ = 0) = 10""^ and 
6 X lO""*, Zobs — 106.6 and 126.2, respectively. In both 
cases, Zobs — 0.3Agw. In the early phase, their ampli- 
tudes gradually increase as their orbital radii decrease. 
The amplitudes reach their maxima when the neutron 
stars merge. Since the binaries change their configura- 
tion to nearly axisymmetric ones soon after the merg- 
ers, the amplitudes drop quickly, leaving quasi-periodic 
waves of small amplitude which result from small non- 
axisymmetric perturbations. Since the equation of state 
is not stiff enough to allow for the bar-mode instabil- 
ity to set in GO], the amplitudes will monotonically fall 
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in the emission timescale of gravitational waves. For 
Pmax(^ = 0) = 10^'^, the amphtude rises and faUs more 
sharply than for /9max(i = 0) = 6 x 10"'', because a higher 
density merged object is formed more rapidly due to the 
stronger relativistic gravity. However, the overall wave- 
forms of the two cases closely resemble each other and 
no significant difference is found besides the difference of 
their amplitude and frequency. This indicates that the 
waveforms only weakly depends on the compaction of 
original neutron stars if black holes are not formed after 
the merger. 

It is interesting to consider the evolution after the for- 
mation of a newly formed, differentially rotating mas- 
sive neutron star of J/Mg 1 such as those obtained 
in the above simulations. Since the mass of such neu- 
tron stars {Adg ^ 2.5 — 3Mq) is much larger than the 
maximum allowed mass of the original neutron stars of 
zero angular momentum {Alg ^ 1.5Mq), the new star is 
strongly supported by the rapid rotation. Hence, if the 
angular momentum is slightly dissipated or transported 
outward, they could collapse to be black holes eventually. 
As discussed in , there are many processes which con- 
tribute to the angular momentum dissipation and the an- 
gular momentum redistribution; e.g., neutrino emission, 
magnetic radiation, viscous dissipation, and gravitational 
radiation. Since each process can affect others in com- 
plicated manner, it is difficult to give a precise analysis. 
Here, we present a rough and conservative upper limit of 
the timescale for collapse to be a black hole. 

Since it will be cooled mainly by the neutrino emission 
|^l| , the new neutron star will first contract on a neutrino 
emission timescale ^ 10 — 100 sec after the formation. 
According to however, the contraction will not lead 
to a black hole because J/Mg will remain ^ 1 even af- 
ter a large amount of neutrino emission plf . Moreover, 
even in the case J/Mg < 1, the maximum allowed mass 
of neutron stars after the cooling may not change sig- 
nificantly if a realistic equation of state is taken into ac- 
count . The bulk viscosity could play a dominant role 
when the merged object has high temperature |43] and 



high non-axisymmetry. Although the temperature could 
be sufficiently high just after its formation, the merged 
object seems to remain nearly axisymmetric during the 
early stage as suggested by the simulations in this paper. 
Therefore, the neutron star could not collapse in t^. 

Nevertheless, the newly formed neutron star is only 
in a quasi-equilibrium state and is continuously driven 
away from the equilibrium state. As a result, the neu- 
tron star will eventually collapse to a black hole either 
because of magnetic fields or because of viscosity. An im- 
portant role in fact could be played by magnetic fields if 
the neutron stars before the merger had magnetic fields 
of ~ 10^^ Gauss. It seems likely that the strength of mag- 
netic fields of the newly formed neutron star could be- 
come larger than their initial strength > 10^^ Gauss due 
to the compression of matter during the merger and/or 
its differential rotation ||l^. The angular momentum of 
the merged object could then be dissipated by the mag- 
netic dipole radiation rather efficiently. The dissipation 
timescale of the angular momentum is estimated as [[44| 
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where B denotes the typical strength of the magnetic 
field. Thus, the newly formed neutron star could collapse 
to a black hole within a year if B is larger than 10^^ gauss. 

If the magnetic field is not very large, on the other 
hand, the dissipation by the shear viscosity may play a 
dominant role for evolution of the newly formed neutron 
star. If the shear viscosity is present, the differential ro- 
tation is forced into rigid rotation in viscous timescale. 
In this case, the angular momentum around the central 
region is transported to the outer parts and the centrifu- 
gal force around the central region is weakened. Because 
it has a very large mass which is probably the nearly 
maximum allowed mass for the differential rotation law, 
the neutron star might collapse as a result of the outward 
transport of the angular momentum. If the microscopic 
viscosity is dominant, the dissipation timescale can be 
estimated as 
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(4.3) 



where we use rjs — 347p^/'*r~^ as the shear viscosity 
of neutron star matter, and T denotes the temperature 
which could be dropped by neutrino emission from lO^^K 
to 10®~^K in ^ 1 year [Q. It should be noted that even 
if the magnitude of magnetic fields is not very large, they 
could be at the origin of an effective viscosity. This is be- 
cause the differential rotation of the star could cause a 
shearing instability resulting in the change of the mag- 
netic field configuration and outward transport of the 
angular momentum j46| . The effective viscous timescale, 
thus, might be much shorter than that of Eq. ( |4.3| ). 
We also have to consider the effect of gravitational 



wave emission due to the bar-mode |40 4^ or r-mode [ p7[ , 
by which the non-axisymmetric perturbation may grow 
contributing the dissipation of the angular momentum 
or the redistribution of the differential rotation. Such ef- 
fects are likely to be important, in particular after the 
newly formed neutron star contracts due to the neutrino 
emission, because in such stage the ratio of the rotational 
kinetic energy to the binding energy could become large 
0.14) due to the contraction |^l[] and/or the equa- 
tion of state could be stiff (F > 1.8 [Q) enough to allow 
for the bar-mode deformations. Since there are too many 
uncertain aspects such as the initial amplitude of the non- 
axisymmetric perturbation and interaction with viscosity 
p8[ , the dissipation timescales cannot be estimated in a 
simple manner. However, it is likely that such effects 
contribute to the dissipation of the angular momentum 
and make the timescale shorter. We can then conclude 
that the newly formed supramassive neutron star (if it 
could be formed) will eventually collapse to a black hole 
in the timescale ^ 10^^ years as a result of one of the 
above dissipation processes. 

The argument presented here suggests that the 
strength of the magnetic field of the newly formed neu- 
tron star is one of key parameters for the evolution. It 



may well be important to investigate the evolution of the 
magnetic field during the merger incorporating a solver 
of the magnetic field equation. 

V. SUMMARY 

In this paper, we present our first successful results 
of numerical simulations carried out using fully general 
relativistic 3D hydrodynamic code. We have performed 
a wide variety of simulations for test problems toward 
more realistic simulations of coalescing binary neutron 
stars and have confirmed it possible to obtain the solu- 
tions for the test problems fairly accurately. In particu- 
lar, we have illustrated that it is possible to preserve an 
approximate quasi-equilibrium state of a mildly relativis- 
tic binary neutron star as well as to perform simulations 
of merger between two corotating neutron stars to be new 
neutron stars stably for a couple of orbital periods in our 
numerical code. These results indicate that numerical 
simulations of binary neutron stars for a long time from 
their ISCO to mergers are feasible. 

Since we could not take a sufficiently large number of 
grid points to impose the outer boundaries in the wave 
zone as well as to resolve each neutron stars precisely, 
numerical errors are unavoidably included in the results. 
In particular, we feel that radiation reaction effect due 
to gravitational wave emission could not be taken into 
account precisely. For a more accurate computation of 
the radiation reaction effect, we will have to adopt a nu- 
merical technique such as a nested grid technique which 
effectively makes the computational region bigger or to 
impose sophisticated boundary conditions at the outer 
boundaries ]2l| ] , unless computational power is improved 
soon. These issues should be pursued in future works. 
We emphasize, however, that besides the error induced 
by the incomplete treatment of the outgoing boundary 
conditions as well as the numerical dissipation due to the 
restricted resolution, the simulation can be performed 
stably and fairly accurately. If we can accept such small 
error (for example, say ^ 5% error in the angular mo- 
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mentum) , the code can be used for investigation of many 
interesting problems even at present. 

In this paper, we have paid attention only to test prob- 
lems. Since we consider binary neutron stars of mildly 
large compaction parameter and small T, the final prod- 
ucts of the merger are not black holes. A black hole may 
be more easily formed for larger F and preliminary sim- 
ulations indicate that this is the case. One of the most 
interesting and important issues in numerical relativity 
is to clarify the criterion for formation of black holes. In 
a forthcoming paper , we will perform simulations of 
corotating binary neutron stars of a large F 2, in which 
a black hole could be formed more easily. 

It is well known, however, that the corotating veloc- 
ity field is not realistic for binary neutron stars because 
the viscosity of the neutron stars is not large enough to 
achieve the corotation |^^. Instead, the irrotational ve- 
locity field is considered to be more realistic one [ pO[ . 
For making a realistic model of the final phase of binary 
neutron stars, it is necessary to perform simulations for 
irrotational binary neutron stars. 

Fortunately, several groups have just recently devel- 
oped numerical methods for obtaining an approximate 
quasi-equilibrium state of irrotational binary neutron 
stars in a conformal fiatness approximation ||^-^, pro- 
viding more realistic models of binary neutron stars. 
Namely, an initial condition has already been prepared. 
In the next step [|9j, we will also carry out simulations 
adopting the approximate quasi-equilibrium states of ir- 
rotational binary neutron stars as the initial condition. 

ACKNOWLEDGMENTS 

The author thanks T. Baumgarte, Y. Kojima, T. 
Nakamura, K. Oohara, L. RezzoUa, M. Sasaki, and S. 
Shapiro for helpful conversations and discussions. He also 
thanks T. Baumgarte and L. RezzoUa for reading this 
manuscript and providing valuable comments. Numeri- 
cal computations were performed on the FACOM VPP 
300R and VX/4R machines in the data processing center 



of the National Astronomical Observatory of Japan. This 
work was supported by a Grant-in- Aid (Nos. 08NP0801 
and 09740336) of the Japanese Ministry of Education, 
Science, Sports and Culture, and JSPS Fellowships for 
Research Abroad. 

APPENDIX A: NUMERICAL METHOD FOR 
SOLVING HYDRODYNAMIC EQUATIONS 

To compute the advection in the hydrodynamic equa- 
tions, in this paper, we use the second order scheme by 
van Leer llsitl as Oohara and Nakamura used in their 
works 



14 5^ or as we used in a semi relativistic sim- 
ulation |53|. As we have shown in Sec. Ill and IV, it 
is possible to perform simulations of many interesting 
problems stably and fairly accurately for several dynam- 
ical timescales using this scheme. If we are interested in 
problems in which shocks play a very important role, we 
should use one of the modern shock-capturing schemes 
as adopted in However, in coalescences of binary 

neutron stars, it seems that shocks do not play a crucially 
important role and that they do not change the results 
drastically We expect that the method employed 

here is suitable for our purpose, although there is still 
room for further improvement. 

Since we solve the entropy equation instead of the en- 
ergy equation, it is impossible to capture shocks without 
adding artificial viscosity. Hence, we add an artificial vis- 
cosity using a method similar to that adopted by Hawley, 
Smarr and Wilson |^ . The following artificial viscosity 
was included only in the case when we studied merger of 
binary neutron stars. 

For simplicity, we include the artificial viscosity like a 
bulk viscosity and it consequently appears in the evolu- 
tion equation where the pressure does. As the form of 
the viscous pressure, we choose 
=r 

-{Svf for Sv < 0, 







(^g60)I 



(Al) 



for Sv > 0, 



where 6v = 25xdkV^ and Cvis is a constant which we 
phenomenologically set below. 
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When including the artificial viscosity, we change the 
pressure gradient term dkP in the right-hand side of Eq. 
(2.10) into dk{P + Pvis)- If we simply regard Pvis as an 



If we write the line element as 



additional pressure, we should add the following term is 
the right-hand side of Eq. ( ^.11 ) : 



r 



(A2) 



However, as pointed out by Hawley et al. the first 
term, which includes a time derivative, could cause a 
numerical instability. Hence, we neglect the first term, 
and only include the second one. 

Cvis is determined in a phenomenological manner. In 
merger of binary neutron stars, two neutron stars of 
nearly equal mass finally collide and a shock will be pro- 
duced. The situation may be crudely compared with the 
ID wall shock problem in which we consider the evolu- 
tion of a fluid that initially has an uniform density and 
pressure, but a velocity field as 

\ f'"tn (A3) 
— Vq for X > 0, ^ ' 

where Vq is a positive constant, and — — 0. The 
solution of the ID wall shock problem in special relativ- 
ity can be analytically calculated. We have performed 
test simulations of this problem, and the value of Cvis 
which well reproduces the analytical solution was used 
in numerical computations. 

To mimic the collision of two neutron stars, we set 
initial conditions as Vb = 0.1 — 0.4 and P/p ~ 0.1 for 
the test. In Fig. 31, we show the results for Vq = 0.4 
and Cvis = 10 at i = 0.4. In this test, we set p = 1 
and 6x = 0.005. The solid lines denote the analytical 
solution, and the open circles are the numerical results. 
Although well known spurious oscillation of p, P, and 

is generated , the numerical results agree with the 
analytical solution fairly well. 



APPENDIX B: RADIAL OSCILLATION PERIOD 
OF SPHERICAL STARS 



We estimate an approximate radial oscillation period 
of spherical stars using the method presented by Chan- 



(Bl) 



then the angular frequency a of the radial oscillation for 
a polytropic star of radius R is written as 
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where ^ denotes the r-component of the Lagrangian dis- 
placement. If we substitute an eigenfunction ^ for the 
eigenvalue equation, we can obtain a correct a. Here, we 
evaluate a approximately by substituting a trial eigen- 
function for ^ without solving an eigenvalue equation for 
^. Following Chandrasekhar |28|, we substitute the ap- 
proximate eigenfunctions as 

C = re"-, (B3) 

— 1/2 

where 6 is a constant. In Fig. 32, we show apc as 
a function of pc for polytropic star of (F, K) = (5/3, 10) 
and (2,200/7r) for 6 = (the solid line), 1/2 (the dashed 
line) and —1/2 (the long dashed line). We note that 
beyond a critical density, becomes negative in each 
model, implying that the star becomes unstable for some 
oscillation modes. Even though we use the trial function 
for ^, we can approximately find the real critical points of 
stability {pc - 1.85 x 10"^ foj. p = 5/3 ^nd - 5.0 x 10"^ 
for F — 2). We also find that the oscillation frequency 
in the Newtonian limit {pc 0) agrees well with that 

— 1/2 

obtained in the Newtonian theory {crpc = 1.377 for 
F = 5/3, and 2.187 for F = 2 ||]). Thus, we can rec- 
ognize that this approximate method is fairly reliable for 
estimating a approximately. 

1 /2 

From the figure, we can judge a ~ 0.60pc for 
(F,pe) = (5/3,10-3) and a ^ O.Olp^" for (F,pe) = 
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(2,3 X 10""^), and the resulting oscillation periods are 
~ lO.Spc and 6.9pc respectively. These are in 
good agreement with the numerical results in Sec. III.B 
and III.C. 
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Table II. The list of the maximum density, total rest 
mass A/*, gravitational mass Mg, J/Mg, coordinate sep- 
aration a, orbital period P, and ratio of the emission 
timescale of gravitational waves to the orbital period Rr 
for corotating binary neutron stars in approximate quasi- 
equilibrium states. 







Mg 




a 


P 


Rr 


6 X 10^"* 


2.84 


2.72 


1.20 


38.3 


902 


9.2 


10-^ 


3.07 


2.92 


1.11 


32.9 


694 


5.3 



Table I. The list of the central density pc, total rest 
mass M*, gravitational mass Mg, J/Mg, and rotation 
period P for rotating neutron stars of T = 5/3 and 2 at 
mass-shedding limits. 



Pc 




Mg 


J/M-^ 
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9.94 X 10-'' 


1.67 


1.60 


0.427 


397 


5/3 


2.77 X 10-^ 


1.58 


1.45 


0.598 


163 


2 
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FIG. 1. Schematic picture for a numerical scheme of time 
evolution of variables from n-th to (n + l)-th time step. 
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FIG. 2. Rest mass M, and circumferential radius Ji as a 
function of central density pc for spherical polytropic stars 
of = 10 and T = 5/3 (solid lines) and oi K = 200/7r 
and r = 2 (dotted lines). The filled circles denote the equi- 
librium stars which are used in test simulations (the tests 2 
and/or 3), and the open circles denote the critical configu- 
ration of stability against gravitational collapse. The cross 
in the figure for M* — pc denotes the relation between M, /2 
and the maximum density for corotating binary neutron stars 
in approximate quasi-equilibrium states which are adopted as 
initial conditions in Sec. III.E and IV. 




FIG. 4. a along the a;-axis at selected times {t = 12.0, 17.2, 
and 20.4) for a spherical dust collapse with zero shift gauge 
condition for the ID results (the solid lines) and 3D results 
(the filled circles). 




FIG. 5. The same as Fig. 3, but with the MD (or AMD) 
gauge condition. The dotted and solid lines denote the ID 
and 3D results, respectively. 



FIG. 3. a and p* at r = as a function of time for a 
spherical dust collapse with zero shift gauge condition. The 
dotted lines are the ID results, and the solid lines are the 3D 
results. 
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FIG. 6. (f) as a, function of time for a spherical dust collapse 
in the MD (or AMD) gauge condition. The dotted and solid 
lines denote the ID and 3D results, respectively. 
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FIG. 7. p and a at r = as a function of time (tpy^) in 
3D numerical evolution of a stable spherical polytropic star 
of pc = 10"^, Js: = 10 and r = 5/3. The solid and dotted 
lines denote the results in zero shift gauge condition, and the 
dashed line in the AMD gauge condition. In the simulation 
shown with the dotted line, wo reduce 7^ by a factor 0.5% 
initially, and in other cases, we give equilibrium states without 
any perturbation. 



FIG. 8. p and a at r = as a function of time {tpl^^) in 3D 
numerical evolution of an unstable spherical polytropic star 
of Pc = 2.4 X 10-^ if = 10 and r = 5/3. The sohd and 
dashed lines denote the results in the zero shift gauge and 
AMD gauge conditions, respectively. In these simulations, we 
reduce K hy a factor 0.1% or 0.5% initially. 




FIG. 9. Xjjns and a at r = as a function of time (tpl^^) for 
a perturbed spherical star of pc = 10^^, K = 10 and F = 5/3 
with a quadrupole perturbation of + mode. The solid, dotted 
and dashed lines denote Krms, 2/rms, and ^rms, respectively. 



26 




FIG. 10. The same as Fig. 9, but for x mode perturbation 
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FIG. 11. The same as Fig. 9, but for pc = 3 x 10"^ F = 2, 
and K = 200/7r. 



FIG. 12. h+ and ftx as a function of a retarded time for a 
perturbed spherical star of pc = 10~^, K = V) and F = 5/3 

with a quadrupolc perturbation of + mode. The solid and 
dashed lines denote those extracted at ^obs = 24.5 and 19.5, 
respectively. 
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FIG. 13. The same as Fig. 12, but for x mode perturbar 
tion. 
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FIG. 14. The same as Fig. 12, but for Pc = 3 x 10"^, F = 2, 
and K = 200/7r. 




FIG. 15. The gravitational mass as a function of p at r = 
(pc) for rotating stars of F = 5/3 and K = 10, and F = 2 
and K — 200/7r. The solid lines denote rotating stars at 
mass-shedding limits obtained from exact equations. The 
open and filled circles denote rotating stars at mass-shedding 
limits obtained using a conformal flatness approximation. 
The rotating stars we use in this paper are marked with the 
filled circles. 
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t=O.OOP 




FIG. 16. Snapshots of the density contour lines for p* and 
the velocity flow for {v^,v'") in the equatorial plane (left) 
and in the y = plane (right) for a rotating star at the 
mass-shedding limit and F = 2 at selected times. The contour 
lines are drawn for p«/p, c ~ 10~°'^^ , where p* c = 0.0122 de- 
notes p* at r = and t = 0, for j = 0, 1, 2, • • • , 10. Vectors 
indicate the local velocity field and the maximum length de- 
notes ~ 0.26c. 
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FIG. 17. a and p(= p,e"'^'*/w) at r = as a function of FIG. 19. The same as Fig. 17, but for J/ Jo as a function 
t/P for r = 5/3 and 2. The solid and dotted lines denote of t/P. 
the results for initial conditions in which no perturbation is 
added and the pressure is depleted, respectively. 
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FIG. 18. The same as Fig. 17, but for Xrms and arms as a 
function of t/P. 
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FIG. 20. Snapshots of the density contour hncs for p* 
and the velocity flow for {v^,v^) in the equatorial plane 
for a corotating binary neutron star of F = 5/3 in nearly 
quasi-equilibrium states. The contour lines are drawn for 
max — 10-°-3^ where max = 0.00305 denotes the 
maximum value of p, at t = 0, for j — 0, 1, 2, ■ ■ ■ , 10. Vec- 
tors indicate the local velocity field and the maximum length 
denotes ~ 0.23c. 
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FIG. 21. xl^s as a function of t/P of a corotating binary 
neutron star in an approximate quasi-equilibrium state. The 
solid, dotted and dashed lines denote Xrms, j/rms, and Zrms- 
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FIG. 23. h+ and ftx 
of a corotating binary 

quasi-cquilibrium state, 
those extracted at Zobs = 



as a function of a retarded time 
neutron star in an approximate 

The solid and dotted hnos denote 
106.6 and 85.3 for = 116, respec- 



tively, and the dashed lines those extracted at Zoha = 69.5 for 
N = 76. 




FIG. 22. J/ Jo as a function of t/P of a corotating binary 
neutron star in an approximate quasi-equilibrium state. 
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FIG. 24. Snapshots of the density contour hnes for p, 
and the velocity flow for {v^,v^) in the equatorial plane 
for a corotating binary neutron star of F = 5/3 and 
Pmax(t = 0) = 10~^ for a merging case. The contour lines are 
drawn for p,/p« max = 10"°'^^ , where p, max ~ 0.00305, for 
ji = 0, 1, 2, • • • , 10. Vectors indicate the local velocity field and 
the length is shown in normalization of 0.3c. At t = 1.62P, 
p* max — 0.011 and pmax — 0.0014, respectively. 
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FIG. 25. The same 

as Fig. 24, but for pmax(i = 0) = 6 x 10~*. The contour lines 
are drawn for p*/p* max = lO""'^-' , where p* max = 0.00143, 
for j = 0,1,2, •••,10. At t = 1.59P, p« ~ 0.0019 and 
Pmax ~ 6 X 10~*, respectively. 
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FIG. 26. The density contour lines for p, in the y = 
plane after merger (t = 1.62P) of a corotating binary neutron 
star of r = 5/3 and Pmax(i = 0) = 10~^. The contour lines 

arc drawn for p,/p, max = 10^""^^, where p, max = 0.00305 
denotes p* at r = and t = 0, for ji = 0, 1, 2, • • • , 10. 
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FIG. 28. The angular velocity along the a;-axis (the 
solid line) and y-axis (the dotted line) at t = 1.62P for 

Pmax(t = 0) = 10^'^ and along the a;-axis (the dashed 
line) and y-axis (the dotted-dashed line) at t = 1.59P for 
Pmax(< = 0) = 6 X IQ-^. 
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FIG. 27. The same 

Pmax(i = 0) = 6 X IQ-* 

are drawn for p«/p« 
for J = 0, 1, 2, • • • , 10 



Eis Fig. 26, but for 
at t — 1.59P. The contour lines 
= lO-''-^^-'', where p. max = 0.00143, 
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FIG. 29. Fraction of the rest mass inside a coordinate ra- 
dius r as a function of tjV for pmax(t = 0) = 10~^ (the solid 
lines) and 6 x 10~* (the dashed lines). 
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FIG. 30. /i+ and ftx as a function of retarded time 
in the merger of corotating binary neutron stars of 
p^^^{t = 0) = 10"^ (the solid hnes) and 6 x 10"" (dashed 
lines). 



FIG. 32. The angular frequency {<Jpc ) of the funda- 
mental radial oscillation as a function of the central den- 
sity for spherical polytropic stars of {K,T) = (10,5/3) and 
{K, r) = (200/77, 2) with the trial functions of the Lagrangian 
displacement. 
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FIG. 31. Numerical results of p, P, and for the ID wall 
shock problem in special relativity. The solid lines denote the 
exact solutions, and the open circles are numerical results. 
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